{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Construction\n",
    "\n",
    "$$\n",
    "\\newcommand{\\sumN}{\\sum_{n = 1}^N} \n",
    "\\newcommand{\\sumn}{\\sum_n} \n",
    "\\newcommand{\\bx}{\\mathbf{x}} \n",
    "\\newcommand{\\bbeta}{\\boldsymbol{\\beta}} \n",
    "\\newcommand{\\btheta}{\\boldsymbol{\\theta}} \n",
    "\\newcommand{\\bbetahat}{\\boldsymbol{\\hat{\\beta}}} \n",
    "\\newcommand{\\bthetahat}{\\boldsymbol{\\hat{\\theta}}} \n",
    "\\newcommand{\\dadb}[2]{\\frac{\\partial #1}{\\partial #2}} \n",
    "\\newcommand{\\by}{\\mathbf{y}} \n",
    "\\newcommand{\\bX}{\\mathbf{X}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section demonstrates how to construct a linear regression model using only `numpy`. To do this, we generate a class named `LinearRegression`. We use this class to train the model and make future predictions. \n",
    "\n",
    "The first method in the `LinearRegression` class is `fit()`, which takes care of estimating the $\\bbeta$ parameters. This simply consists of calculating \n",
    "\n",
    "$$\n",
    "\\bbetahat = \\left(\\bX^\\top \\bX\\right)^{-1}\\bX^\\top \\by\n",
    "$$\n",
    "\n",
    "The `fit` method also makes in-sample predictions with $\\hat{\\by} = \\bX \\bbetahat$ and calculates the training loss with \n",
    "\n",
    "$$\n",
    "\\mathcal{L}(\\bbetahat) = \\frac{1}{2}\\sumN \\left(y_n - \\hat{y}_n \\right)^2.\n",
    "$$\n",
    "\n",
    "The second method is `predict()`, which forms out-of-sample predictions. Given a test set of predictors $\\bX'$, we can form fitted values with $\\hat{\\by}' = \\bX' \\bbetahat$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LinearRegression:\n",
    "\n",
    "    def fit(self, X, y, intercept = False):\n",
    "\n",
    "        # record data and dimensions\n",
    "        if intercept == False: # add intercept (if not already included)\n",
    "            ones = np.ones(len(X)).reshape(len(X), 1) # column of ones \n",
    "            X = np.concatenate((ones, X), axis = 1)\n",
    "        self.X = np.array(X)\n",
    "        self.y = np.array(y)\n",
    "        self.N, self.D = self.X.shape\n",
    "        \n",
    "        # estimate parameters\n",
    "        XtX = np.dot(self.X.T, self.X)\n",
    "        XtX_inverse = np.linalg.inv(XtX)\n",
    "        Xty = np.dot(self.X.T, self.y)\n",
    "        self.beta_hats = np.dot(XtX_inverse, Xty)\n",
    "        \n",
    "        # make in-sample predictions\n",
    "        self.y_hat = np.dot(self.X, self.beta_hats)\n",
    "        \n",
    "        # calculate loss\n",
    "        self.L = .5*np.sum((self.y - self.y_hat)**2)\n",
    "        \n",
    "    def predict(self, X_test, intercept = True):\n",
    "        \n",
    "        # form predictions\n",
    "        self.y_test_hat = np.dot(X_test, self.beta_hats)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try out our `LinearRegression` class with some data. Here we use the {doc}`Boston housing </content/appendix/data>` dataset from `sklearn.datasets`. The target variable in this dataset is median neighborhood home value. The predictors are all continuous and represent factors possibly related to the median home value, such as average rooms per house. Hit \"Click to show\" to see the code that loads this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "boston = datasets.load_boston()\n",
    "X = boston['data']\n",
    "y = boston['target']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the class built and the data loaded, we are ready to run our regression model. This is as simple as instantiating the model and applying `fit()`, as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = LinearRegression() # instantiate model\n",
    "model.fit(X, y, intercept = False) # fit model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's then see how well our fitted values model the true target values. The closer the points lie to the 45-degree line, the more accurate the fit. The model seems to do reasonably well; our predictions definitely follow the true values quite well, although we would like the fit to be a bit tighter. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "Note the handful of observations with $y = 50$ exactly. This is due to censorship in the data collection process. It appears neighborhoods with average home values above \\$50,000 were assigned a value of 50 even.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEqCAYAAAD6aUxzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deXxU5b3wv89sySQBEkJYFFCkiKYYhCiy3NuiqNWWai2orYJrBaQWr6916b2X2vdSe0X02qpVwLYq4gKCXn3popZKF5SqEbdGERcsKJAQEsgyme087x8z5zDLmayTmUny+34+fDJzzplznnOG+f2e57cqrTWCIAiCAODI9gAEQRCE3EGUgiAIgmAhSkEQBEGwEKUgCIIgWIhSEARBECxEKQiCIAgWohQEQRAEC1EKgtDLUUrdoJS6IdvjEPoGrmwPQBCErqOUWgj8d/R1s9Z6dZaHJPRylGQ0C0LvRCl1HPA28H+IrPrvAiq01p9mdWBCr0aUgiD0QpRSDmAL8LHW+srotjXAMcDpWmsji8MTejGiFARBEAQLcTQLgiAIFqIUBEEQBAtRCoIgCIKFKAWhz6KUulEppZVSN6bYP14p5VdK/aWNc0yLnuOZNo55P3qewTHbzlNKbVZK7Y3u+0Ip9Wel1OJs35MgtIUoBaEv87fo36kp9t8HOIHrUp1Aa/0qsAOYrZQqTdyvlJoCnAD8P631wei2BcBzQDnw/4C7gd8BXuDKLt3JEbp9T4LQFpK8JvRl3gR8wGmJO5RSFwJnAfdqrd9p5zyPAj8Dvgvcn7Dv8phjTBYCAWCi1rom4bpDOjx6e9J1T4Jgi4SkCn0apdSfga8AR2utv4huKwQ+ADzA8VrrQ+2cYyTwGfCm1vrUmO0eYC8Qip4/FN1eRWT1MFJrXZ+L9yQIqRDzkdDX2Rr9G2tu+TEwErilI8JTa70H2AycopQqj9n1TWAw8LipEKI8DhQA/1BK3aOU+pZSqqw7N5FAt+9JEFIhSkHo65gC9DQApdQJwA3Aq8SbfNrjkejfy2O22ZmO0Fr/T3TfP4ElwLPAfqXUy0qpUzoz+BSk654EIQkxHwl9GqVUCVAH/FVr/VWl1B+B04FTtNbbO3EeL7APaAJGE1khfAH8Q2t9chufKwamAxcAVwENwImJvobOkK57EgQ7ZKUg9GmiNv33iZh+LgFmAas6Kzy11j5gPXAUcCZwKZFAjTZn5lrrBq3177TW1xBZbQwG/rWz95FwzrTckyDYIUpB6A/8jYiNfxVwAPjPLp7nkejfy6L/QkT8B3Eopc5RStlF9g2N/m2JOXasUuoEpZS7k2NJ1z0JQhwSkir0B7YCC4Ai4AYzn6CzaK23KqU+Ai4E3ERyE+zMQE8BrUqpvwG7AEVkdXAqUAX8MebYzUQqm46JHttR0nJPgpCIrBSE/oDZX+B14NfdPNejRBSC+dqOW4k4fScDi4kkrLmBW4iUtQ52cwyQ3nsSBAtxNAt9HqXU88A3gKla69ezPZ500BfvScgNZKUg9GmijthvAg/2FeHZF+9JyB1kpSD0OZRSo4FLgLFEHMI7gSla65Y2P5jD9MV7EnITcTQLfZFziDSzbyBSmO7f+oDw7Iv3JOQgslIQBEEQLMSnIAiCIFiIUhAEQRAsRCkIgiAIFqIUBEEQBIteH310zjnn6D/84Q/ZHoYgCEJvQ9lt7PUrhQMHDmR7CIIgCH2GXq8UBEEQhPQhSkEQBEGwEKUgCIIgWIhSEARBECxEKQiCIAgWvT4kVRAEoT9hGJq65gCBUBiPy0lpoQeHwza6tEuIUhAEQeglGIZmx/5GrlnzBnvqfYws8fLQZacwftiAtCkGMR8JgiD0EuqaA5ZCANhT7+OaNW9Q1xxI2zVEKQiCIPQSAqGwpRBM9tT7CITCabuGKAVBEIRegsflZGSJN27byBIvHpczbdcQpSAIgtBLKC308NBlp1iKwfQplBZ60nYNcTQLgiD0EhwOxZeGFLJuwVRChsblUAwtypPoI0EQhP5IKGSwo6aJRWurrOijlfMqOWHYAFyu9Bh+xHwkCILQS6hp8lsKASJO5kVrq6hp8qftGqIUBEEQegnBsGEbfRQKG2m7higFQRCEXoLb6bCNPnI50yfKRSkIgiDkAIahqW3083l9C7WNfgxDJx0ztCiPlfMq46KPVs6rZGhRXtrGIY5mQRCELNPR8hUul4PxQ4uSoo/S5WQGWSkIgiBknY6WrzAMzUcHmrl49Ta+umILF6/exkcHmm1XFV1FlIIgCEKW6Wj5ikzUPhLzkSAIQpYxy1fsqfcxaVQxi2aOpbTQg1IKw9CWCSkTtY9EKQiCIGQZs3zFPS/t4PLpY7hl4zu2vgWPy8nCfz2WuaeMxulQhA3Nhjf+mdbaR0rr9NmissEpp5yi33jjjWwPQxAEoVsYhmbf4VYuWvVq3GpgZImXZxfPoGxAHsFgmA9qmrg2JqP5wXmVnDC0CLe704rBtjaG+BQEQRC6SUfCSdvD4VBords0D9U2ByyFYO67dm0VtX2ln4JSyqmU2q6U2hR9P0Yp9Xel1E6l1DqlVPpK/wmCIPQAZjjpBQ9sZcbyl7ngga3s2N/YJcXQXmnsVBnNwT6U0Xw98H7M++XAPVrrcUA9cHVWRiUIgtBB0hkR1F5pbJdD2Wc0p7FKataUglJqJPAN4FfR9wo4A9gQPeRR4FvZGZ0gCELHSGdEkMOhGD9sAM8unsHWW07n2cUz4hLYivKdPJiQ0fzgvEqK8tPnaM5m9NHPgZuBAdH3pUCD1joUfb8HONrug0qpBcACgNGjR/fwMAVB6MsYhqauOUAgFMbjclJa6OlUf4LYcFKT7nRDczgUZQPsy1b4Q5r7Nn/I0tnlFHvdNPiC3Lf5Q356wUldupYdWVEKSqnZQI3WukopNdPcbHOorVFOa70aWA2R6KMeGaQgCH2ejpaXaAvT5JN4jnR2QzMJhgxerK7hxeqauO23fTN9PoVsrRRmAOcppb4O5AMDiawcipVSruhqYSTwRZbGJwhCPyCVP8AMAU1F4upiXFkRzy6e0eXVRkdJ96rEjqz4FLTWP9Jaj9RaHwt8B/iT1vpS4GVgbvSwy4HnsjE+QRD6B13xB9hFG+2sbaK00MPRJQWUDUhve8xY+mOP5luAp5RSPwW2A7/O8ngEQejDdHTmHbsyUEp1aXWRDmId0T21Ksm6UtBabwG2RF9/AkzJ5ngEQejddMZx3BF/QKLfYcOiaT1ef6gt2nJEp4OsKwVBEIR00VnHcUdm3ol+h7rmQI/b9bNJtpPXBKHLpKO0gJAZMvVddSWRzJx5m/4AIGGs8VnEK7d8zPI5FT1q188mslIQeiXpCCUUMkMmv6vuJpLZjXXV/ErOLh9qhYFu393Ao698yvqF09Ba92i0Uaoxdievoj1kpSD0SjLRbERID5n8rlLVDlJKJa1S7FYvdmNd+FgV//mN8riVwQ1njWf4wPy0RRt1dCWVzjpLqZCVgtAryUSzESE9ZPK7snMcr5xXyU+ef48Xq2usVcq4siJ21jYlrV4GF7htx+p0qB6L+OnMSqqreRWdQVYKQq+kvWqSQs/T0dltJr+rxNpB6xdO497NH1qmH1OI1jT5bYVrWJNyrLF+h3SaazqzksqEghWlIPRKMpHEI6SmM2aMTH9XsY5jrXVSSQiz1LSdcNVaZ/z/VWcEvVL2VVIj9UTTg5iPhF5JJpJ4ejOpnJHpclJ2xoyRze8qVXKa2+lIGVY6fpg3o2PtTOkKp4Llcyri2nUun1OBM43DE6Ug9Fp6Oomnt5LKRp3Kjt6VKKDE2a3ZbL4lEKK2kSRB2tZ31ZPRNKmS04YW5aVMWsv0/6vOFNRzOBz8Zcd+Hr7i1LgezccP+1LaxiM9mgWhj1Hb6OeCB7YmzTzXL5zGT55/jytnjGH4wHzCWnOgKcDxw4oYXNjx4m+lhR7qmgPWNSaNKuaHXxufstl8W2QiXLWnV02dvW53jg2FDD7Y38iimB7NK+dVcsKwAbhcnfYG2A5GlIIg9DE+r29hxvKXk7a/euvp7G/04wuEuWnDEQG+an4lJw4faCuEDEOzq66Zz+paKPA4aQmEOaa0gNElBdaqY+nscpZtqk7ZbL4tUimwTNQR6kl6Stml+XnZDkQczYKQo3Q1CzhVtI+hob45aCkEOBKHnypnoMEXYP/hVpY+9x4Xr97G0ufeY//hVg77g5af4IThA7ocEZNLocXpzLruqdwMiT4ShH5Kd5KUUkX7OB2KUYO9tkLFFwjZCsLYVYV57E0b3sEXCFu29wKPq8shp+2Fq2aqPEa6k8J6SnhnIrxXHM2CkIN0J0nJLtqnxOtmR00jNYf9tpEu7+9rZNmmasshXe8LEgiFCRnaVriFY2RldzqPtfXZTJbHSHdSWE81wykt9LDmqilJ5ry+3E9BEATanml2xCmZGEFT2+hn4WNVlBXlJYU03n3hRO74/QeWIHzie6dxya/+zp56H08tmGoJNzPCqLTQgyvqqHU4VJshp+2Nta3P1jbaJ5j1hL8h3TP7nmzR6Q8ZLH3uvbjzphNRCoKQg6SaaXo9zi7Nnk2ht6fex10v7LAavw8bmM8N695i++4GICIIaxr91nWX//4D7rloIg/99RMunz4mZYSRXRhnR2f6qUJAM+lvSPfMvqdyM6TMhSD0U1L5BUKGbtOBmcoGH2uL3r67gYWPVXHj02/z6YFmSyGY14l1hm7f3cDPfvcBS2d/2VIIdte1I5UA23e4tUM+gkyWx+iJrOvEktzpMHllQlHKSkEQcpBUM829h3xtmpVSzcxTmTPyXEcye83w1F/88cO489c2+fF3QRilEmBfNPiYu/LVdlc5PWmCSaS3ZMj3lK8iFlEKgpBBOpPQZGdWSVm2weVg3+HWNk0LdkIPSHJI33DWeKr3NsaVUTjQ1PluY6nGaq4u2jN9dFZQdzcZrTdkyGdCUUrymiB0g85mrXY3msb2HPNPwe1SNLQEmbvy1aTPbL3ldI4uKejUPe073Epto58Cj5PWYJgmfwiXQ3HD+ret6665agpF+S6CIcP23u3GunxOBXe9sCPOZNXZ8XX4ufTRpktpzMSWjGZBSCedFUTpykaNFQpul4Nw2ODCVdtSZhavXziNo4q93VZgsUrA63Gy/7A/bv+qeZWMKM6n2Hsk8qjBF8AXCBPW4HIoq69Bd+7fjr6aGd3DiFIQhHTSVo2h4QPzk2bNpj+gwRdk5ZaPrdlyWzPl9gT5wWY/h3xBag77MbTGoRQ3Pv123Mz8S2WFlA3I7/RMuq1rp7r3ZedPYPigfNvie2uumoI/ZPTIbD5VaY90rEJyjZ5eKYhPQRC6SFuO1EO+oCXs2jKj1Db5k+zy5o/eMAwONAdY+FiVrRA1DM3ehlYWxhRHu/+SSdzx7ZPIdzsZ5HWz4oUPuP2CipSRQM8sns7QAfm299eWjT3VvRd4nFyz5g3WL5yWdL3LfvMaz183o0ecuZlwwOYCmTCTSUiqIHSRVCGTwbDBvkOt7GmIhF0eaE5Owrpl4zssmTUuyUkYW27hrT2HLIVgfs4MAzXt/qZCMPdf98R2gmFN2NCseOEDbjhrPKWFnpRCvMUf7lIph1T33uALsqfeRyhFExtfINwjHcz6S9OlTPS7lpWCIHQRu0iQ+y+ZhD94JOP07PKh/Mc3ym0F5JeGFpHndrD3kC+uJLV5vmKvfb/gQCjMjv2NNPtDtvtHDy6gKM/J7RdUWDPxVDPpTw80U5jn6rTd3e7ezdXPyBIvrjaa2PQEvSWktLtIQTxByGFMQbR+4TQ2LJrG0tnlNLWGLJv+pFHFXD59DJ/UNtvOqrXWfPuBV6wCbO/vO4ziSK2hBl8wZetFc3Zot7+msZWQoQmEwtaqorTQw6p5lXEz6eVzKrh3884uCRTz3p9ZPJ0tP5zJsvMnWOaw2CY2mZy590SyWK6RiYQ+cTQLQjeJtfPefeFELl69DYBV8ytZtqmasqK8pCY0ZpJYYiTOk9dMZdmmf/BidU3K5jXFBS6m/ffLTBpVzE/O+zLff+LNuJVKMGTEhY6aNucGX4C3dx+iwOO0nN21Tf5uR+hkq4lNfyTNPgWJPhKEniJWAF68eht76n2sWzDVUhBmMblir5uRJV4cDpj238nRMpv/z1dxORWXRgvSmeYnl0OR73HQ4jfwhwx2H2zh3s07OWN8Gd+aPJJg2CBsaJr9Ia57crttaGZpoadbAkWEfG4g0UeCkCG682MzTReGoS1bu2n+2VPvY/vuBlZu+Zgls8YR1hoMxdnlQ5NWCg4FDS1Bq2Bdgy/I7b+t5o45Few91BoXiXTPRRNxuxx896Ft1rbHrp5i7+ANRkxEXbW7d2aGKsqjZ+npzGtRCkK/oS1hla5leazD0zAMVs2vtEpW33zO+Lg2mA9cOhmAF6trGFniZcXcSDmJgy2RMNRYbvtmOCkS6Yb1b7Ps/Alx23YdaLF18H7R4AOtcUadziMGedu8r8RnpbEvxJdoeupPmcV9FVEKQr8gFDLYUdOYMuY/nSWJY2dyZQPyrZm5aVYyz7/48TdZc9UUbj7nBJwOhcuhCBmah/76cdz5Iq007ZvdFHjiHYz3bt7JqnmVVqjq2eVD+dG5J4KCj2ubuXfzTssZnEpQ2wn2tVef1qGol0yUdhZ6lqwoBaVUPvAXIC86hg1a69uUUmOAp4DBwJvAfK11+gJwhX6JYWi+OOSzjfk3hVWqUD9fMGw1k+kKpoL4vL7F9vxKwRW/ed0Svj+/+GQWn/6luIJ0a66agtbYrgBaAvFCubbJz4jifNYtmMrB5gCGhvm/eS0pbLQtQW0K9rKiPMuMFTK0rbkrMeoll3ouC10jWyGpfuAMrfVE4GTgHKXUVGA5cI/WehxQD1ydpfEJfYi65kBc4xiTSDJVpDexO1pCOpaRJV4+rmnqVq9ek1ShhLsOtMQpqn9b9xZNrSGWzi7nLzfN5NnFMyjKd/HT31azfE5FXIjnqvmVHFNakBT2Wez14HE5OdAUsCKTzPPfsvEdFs0c26agDoTCVsTUsk3VXLx6G1c8/Bo/mHU8Z5cPjbtWYohpJnsgdJVM9X3urWRFKegITdG37ug/DZwBbIhufxT4VhaGJ/QRzB9/SyDEIK/bEmgmI0sivYkveGArTa2hpLh6M44/HRmjdhm3D146mXs374w7bk+9D7fTwbJN1RTlRxbyvkCYOZWjeG775yydXc66BVNZOrucIYUeji0t5NnFM9h6y+k8u3iGZRIqLfQwZkihrSI0I6BSCWqPy8mSWeOSmupcu7aKn5w3Iela7d1nLmUWx2aMm/kh6VD6fYms+RSUUk6gCvgS8EvgY6BBax2KHrIHODpLwxN6OXZ28VjH7tnlQ/n3r5dT3xJg6exy7vj9+9wxp4J1C6ZaRetiSzx31/yRmHGrlMKpIuaeWEyT0JqrpiRVIY0tO22GmaaKRHE4FAV59lnMLYFwm4K6LYWitW6zwFyuZxaLz6N9sqYUtNZh4GSlVDHwLHCi3WF2n1VKLQAWAIwePbrHxij0LmIjZpRS3PPSjiTH7hPXRGbZh30h5v3675bAfXBeJcGQAWBlJJuky/yRKMBjw1cTy0+HDM1lD7ySZPoxy2N3ZPY9pDAv5fnN8tapxplKoXTkOeRysxrxebRP1qOPtNYNSqktwFSgWCnliq4WRgJfpPjMamA1RJLXMjVWIfeIVQRhQ/PT31ZbIZ7L51RQ2xiIb0p/uJW65kBc3wHTNLLs/Ancu3knK+ZWxIWO2gngdMTitzWrTuWYPnH4ACsRrb3rdWfWbqdQcskM1FX6SzXV7pCt6KMyIBhVCF7gTCJO5peBuUQikC4HnsvG+ITeQaqS1MVeD7PKh5HncnDn3Apu3vCOZXKpaw6kLDRnhneGDc0jV04hz+Vg36FWhg2Mr6PT1Vj8VIrEbladSnh5PS4rSa620d+usLc7f0cUWq6bgbpKJvs+91ayUuZCKVVBxJHsJOLsXq+1/i+l1HEcCUndDszTWvtTn0nKXPRn7Bq9nF0+lO+fPi6uHtDdF07EoRSDizzUNwcoG5BHbaOfuuYAm6v3M6t8GKWFHsoG5LHvUKttk5phg7xtXretLl9mBzKz98H040pZ8NWxuJ0Kt9PB0KI8XC5H0mdSKR6gywliklwmGdcxSO0joW9h123LLEKXKLCXnT+BKx953cocvvMPOygb4OG6M8ax+PGIAnn4ilOtktexnzXbWbZ1XbDv8mUK4X2HWln63HtMP66UedOOsa45ssTLynmVnDBsgK1isBNe3Wk9KW0rhRhslYKUzhZ6LXYx8aWFHvswzAK39fqmDZFY/TmVoyzhDFDgcaaMuGnvuiNLvLidyT8nM9rFPPc1Xzku7pp76n0sWltFTVPygjhVKejuOEvT7WiVmP++hygFoddiFxNfNiDPVmAX5R1xn00/rpTyEQMZP3wAD19xKjeeOY5V8ysZmuKziU7I0kIPD82Pv+6KuRUcbg0SikYwmZhC2CyO53QoW6EcCsd/ri26kyCWzuQyw9Dsqmvmvc8Psafex3ufH2JXXbMohl5O1qOPBKGrOByKsaUFrFswlZChcTkUXo+D5XMq4noQLJ9TQWu0SuiNZ45j5onD4iqLPjivkvs2f0htYyAp8mjlvEpKvO6k65YWeVh2/gSrN4Fpjrrtm18GsMw9phBeueVjls+pwEhRrsJls8pIRXecpel0tDb4Auw/3GqZ3EzlWFzgZnChmKJ6K+JTEHKW9hyCwWCYD2qauHZtpArpklnjOHZIAY2tIQ77gridDoJhgzyXk5JCD/sO+RgzpDCuMB1EhPLS2eUsfKyKSaOKWTJrHKNLCzAMTSBsUFLgYfjA/Lhrf1bXzFdXbLHep2qIM66siJ21TVYtoZ9+awIGcO3aqnZ9Ct15Nj312Vg+r2+xfZbrFkxtM8FNyBmkn4LQe+hIlExNk59ro9E8i2aO5WBzgA/3N7GxajdXzhjDxqo9XDD5aG5Y/9aRiqFfPzFl6QeA7bsbuPKR1/nzTTO5/ffvWzkPidd2J/QgXjRzbFJZCDNTNjG0szjfxfqF0wiFDVwpoo/aozsJYulKLgunqNwa7t3zzH6P+BSEnCRVOYLYGkQhQ1NWlMelU4/hst+8xtyVr7JsUzWXTx/Dw1s/5QezxlmmIIA5laOsfgOxjCzx0uALxr3/pLaZy6ePYdKoYttrDy3KY2VMz+NUDu5AKJzkMHa7nRxV7GV0aSFHFXs7rRByhXy3vX8i390770eIIN+ekJN0JErG5VAsmTXOthLonMpRST0Iir1u7t28M6na6Mp5lWys2m29NwvhmRVFba/tcnDCsAGsXziNv9w0kxGD8nO+Omi6MbOeE4vfDRF/Qq9GzEdCTtKRcgRDi/LQ6LjWlSu3fMz23Q2UFnoIhnXcORp8QWqb/Nz1wg7rMy2BMAUeJz/+5pe5+l+OSyqEZ5qV7AS8y+Ww8hfs6hj19UzZvpr13N8RpSDkJB2JknE6HRzyhaxkNXOW/+grnzK40MPKLR/z4KWTuTaaF7Cxarf1fuFjVVblVF8gTL7baVsIzwwlbU/A91cBmcvF7/oqPZ2RLdFHQs4SChnUNPkJho24chCxPwq76Je1V5/GAy9/xPqqPTy7eBoNLSErdPTDvYf5xsSjONgcoK45wMaq3SyZdTzHlxXyYW0zixKigsqKPDgcjqQfXk/+MKUMg5CKNJcpkegjwZ5sCKH2rmkY2grlTBXiefeFE239DvUtAdZX7WHSqGIcSpHvdlirgKcXTuOzuhar+F1tY4BFa6tYv3Aa927+kDu+fRLDB+XjVIoDTQHcLkdSzH1P1g+yO3dHyl0L/YNM9IMQR3M/p6c7UdmVQejINVP9569pOtJ4xjTtxDKyxEtDSySSaMmscfzy5Y+iVU9P5c83zaSk0M3S597j4tXbWLapmh9+bTxlRXkEwwa1jQGK8l3sPuijptFPkz/EweYAB5v9HRpbqu5sic8gFDJSloawO/fCtVW8vfuQdAgTMtIPQlYK/ZyenHmkmlEPG5jX7jVT/ecPhg1ru5klHJswdveFE/F6nPz5ppm4nYqr/+W4uKqnv/jOyayYW4FDKRp8QR595VOWzBoHwIoLK2gJhJMydOubgxTl+62VQGd+mHbPYOW8Su7d/GFSDgSALxhKWdZbOoQJmegHISuFfk5PzjxSKRxfoP1rtlV0zty+fXcDd72wg2XnT+BPN36VpbPL+fXfPqHA42TvoVa0ju+itqfex/VPvYVSylopXD59DF8aWsjL7+8j3+0kEDJYOrvcyk+4acM7DCnyxK0E2qsfFLsy2He4NekZLFpbxZzKUXHPpMEXYMf+Rj6uaU6ZR5HuGaHQ+yjxuuPyY1KVYukOnVYKSimvUuojpdTOaIMcoReTzgJpiaRSOOFo/Z+2rpmqAbzLASvmHskzqG3yU1rkYWC+iwlHDeS/55yEP2jww6ffpjVo2F5/SJHHen3LxncAmHRMKd9Zvc1KgPvh18ZbisHM3DUFclvN6RNNY180+NrMoDbftwbCXLPmDds8iuVzKli55eM+n/cgtE+9L8i9mz9k6exy1i2ItJa9d/OH1MckX3aXrpiP/gt4K/r6x8CP0jYaIeP0ZCcqu6Xu2eVDcTkUa68+jU8PNHPv5p3UNvmTrmkX4ul0wCe1zTz75uesuWqKFUF0/592cv2s4xlS5CEU0jT6Q9x94UQ8LoftUtupjjhrI/uUbQKc2Q+5oSVoCWTTQT4wWqrCqYiLTqpt9MetDOqaA7ZjSMygDhoRxbOn3sddL+zgjm+fxIhiL/+sa+GuF3bYPiOh/xEIhXmxuoYXq2vitt/2zSz5FJRSk4A5QGV0U5VS6imt9dtpG5GQUXoyvj5R4ZxdPpQfzDqei1a92m5kjV100t5DPuqaA5x70ggu+81rcYK2em8j9333ZFqDhlXa4uzyoTxw6eS4hjYr5law73Cr9bmRJV78IfsVRWmhhwcunYw/aPDQZadQ4nW3G3WUuDqy83uYPgXz+pEVkLKUx/bdDcz79WucXT6U2775Ze6/ZJKEpgpAZnwKkqcg9Cixwl0D37HJK3hm8XSGDsi3jjUMg6ChqWsKsJmrJXUAACAASURBVO9wKxurdnPDWeMpLfTwn//7LreceyKz7v5z0rX+evPpVklsk7PLh7J09pfxh8IcaAowuNDNXS/ssJy8K+ZWEDY0tz7zbtK4nrxmKvdt3sl1s8YxsjjS37m9rmWpWoT+5LwJaK3xuJyUeN3U+4JxCq/BF2DHvsa4st0r5lYwfvgAKUMtWEiegtDrMTNeaxv9fFbXbAnLSaOKWTRzLMVeN4GgQTAY5qMDzdzz0g4unz4mbmb9y0smY2iNQ8GPvn4iLoe9WSix5IXZfzkYNth90MdRxfm88O4+ls7+Mgu+MpZhA/NZ8uR2gKTZ/C8vmcx9m3fyyid13HTOCR2OOrIzx91w1vik0tuJEUTFXg/DBuZbPRpaAmGGDYysoNKBJMT1HfJcjrj/J3lpLqgoSkHICIFQ2LKvlxXlceu5J8SFiq6aX8kv/vghcypHJZWg/v4Tb3LHt0/isC/ITRveYcXcirjyFSNLvNx/ySQaWoJxJS8euHQy9/9pZ9yqYOrYIQTCBsUFbpwORW2T37LjL51dTmmhh8GFHm5c/3aSHb8jS/eumuMcDsWxpYUMyHenXXD3ZLKdkFnqmgNJptPE1Wp36ZD5SCn1ElCktZ6WsP0k4E3gcq31E2kZUScR81HvoLbRz388+w4/OGMczYEwP4wqBHPFUFrooaTQQ4s/xI+f+4e1ijCL3K24sIIrHn6dPfU+1i2YiqE1rUHDKl9R6HHamoDM5jnm+yevmcqyTf/gxeoa/nbLTPY2tHLD+vg8hpHFXkKGThLMvVW42pm00i1IhMzweX0LM5a/nLR96y2nd6WxUbfMR38D/l0plae19gMopRTwAPBKthSC0HsoLfRww1njaWoNMXxgvqUQEruVPXDpZH5yXjnff2J7nG093+WIq3b65q46Lp5yDAejuQMD8t0dCv00tKa2MfKZ1qDBz373QZzJad1ru/nBrHGEDY3bFT9hSlwFKKVwqsjsLZfNMZnIghUyQyYczR1VClsBDzAJ2BbddhkwFZicttEI7ZIrtuHYcXg9TkKGJhgy4l7Hjs8UqA2+AM3+MCNLvLbdyhY//iaPXDmFuy+ciKF1RDg7HWgUZ5cPpbYxwPCBecw+eaS1jB5Z4uXx753WodDP+uYAi2aOZeFjVew71Eptk99aSUwaVczN54yP69+cuBJwOBSlhZ5etWLIhCARMkNPhpCbdFQpbAPCRJTANqVUMXAncL/W+t20jUZok1wxX8SOo6woj5vPGc9NG96Je203PodDMbgwD00rK+ZW4HY6bGewDS0B7vj9B9x8znjLJGT6HQYXevCHDD6tbaasKM+K7b/9t9WsnFcZV+X0wXmV3BcT+rlibgUoOGpQPgB3v/gh91w00TIfLUno1GZmG5vRUSZmpnZZUZ61yth3qJVhA/NyMlIoE4JEyAyZKNHeIaWgtW5SSr1NRCkA3A4YwG1pG4nQLpmokNgW5urAFwyx71ArZUV5LJo51hKkS2eX2wrVxPENyvNE3yvbGWxddDafeK6Fj1Wx7PwJXPnI61aE0OPbPmN91R5erK7hB2eM464LJzJsYB67DrRgGAbfnXKM1Tznzj9EksAeuXIKYGZD5/H0wml83uBjSFTJxLKn3kdr0IjbFgiFKSvKSzJ9rZpXmZOVTPtrr4e+Sk/3sOhM9NFW4Dyl1GRgERHn8uGeGZZgRzZtw3arlOVzKhiY77LGVOy1t+ub44s1Od35hw+4csaYpNn98jkV3PXCDm4994SUheHM199/4k3WXDWFnTVN1Db5+eJQJCnNdGKvWzCVKx95Pele8lyKrbecbglHiPgpNNpWSTkTZKfH5WTJrHFJpq+Fa6ty1nkrzXCEjtKZANe/AccAa4CtWuu1PTMkIRU9WacoFWZxtz0NLdbqAI6Ugoht3p6qlLXH5SQUMthT38Jndc0EQpEy1d996O+8/P5+Hv/eaWxYNI2HrziVR1/5lO27G1KXxY7xEeyp93GwOcCSWeOs+kCxiinVObSGEYO8lA3Ii/N3FHpccXWVTJOT1xP/fEsLPYwZUijOW6FP0tmVAsAJiHM5K2TaNpxqdWD2MN5T76PJH2LF3Apu2vBOJHQ0+jp2fCVeNztqGln4WIy9/9JKBnpduBzK6p5mRiNV721k5ZaP4+z9ZrhoIGSwbsFUGnxBNlbtZpDXHTX7tABHFMGeep9tiYnlcyr46W+ruf2CCmvmbK5ggmGDY0sLuOvCiShImUDmcCgK8sR5K/RNOlzmQik1CKgBHtRa/1uPjqoT9Lc8hUxGH6WKbzdj/0eWeFm3YGq70UepykMsO38CQ4o8/Pi5f7B9dwNwJNO54uiBNAfC7D4YMRlpoLTIw5XRXAXTkbzprT2s+usuS+D/Zcd+Zk8cybWPV1n1j24990QO+YLUNPpZueVjtu9usOK6bTudza9kSKF9G87Y7yEXnP6C0A26Xebix8BBxLmcVTJpG07lwyj2ui0hOGKQt10hmOo8BR4n1z7+puU8hkiPhGWbqnnqmqlWshrAqvmVlq/A/Py1a6t4+IpTOePE4VbDnB9/88t4PQ7WL5yG1hqlFD95/r24qpKxM/q65gD3vLQjLlfhF3/8MG4lYYc4b4W+SptKQSlVAEwE/hW4HrhQa30oEwMTsk+q0tfHlhbwVDSr+ECznyGFeSln02ZjmlQ5BHvqfYwuLbD2mwlssR3WILUT+5AvyMWrt1krhabWIAo3IwZFopj8oTD/8Y1ygLhOZyVeN7WNfoLhMN8/fZxVOts8j2HERxzZIc5boS/S3krhTOA54HPgeq31sz0/JKE7GEZEULcGwziVwutxMjAvuSpnW0LcMAzCGpwqMkM3fQFnlw/lh18bzycHmuP9BvNPYfzwSDtJ07Tldjloag1x2W9eY/pxpUklrE3fhNnf4IlrTsOpFCFD44hmCscqklhfgYkZvgpHHN+PXDkFr8eZZNr55SWTuf7M4xmQ52LEwHx21jZxzZo3rJ4JsSuQWza+w/qF0xCE/kibSkFr/Twp7E7dQSk1ikgU03Ai+Q6rtda/UEoNBtYBxwK7gIu01vXpvn5fxc7Ofc9FEykp9FimmFS2b/OziVVKn144zarIOLjQE8lHiPYwhmguwmNv8PSiaTS0BOOuvWJuBWVFecwqH8b9f9pp2zhmxdwKHnj5Iy6YfHScovnlJZPiwlU3Vu1OKoJnKhaTPfU+3M6IYknM5zCL6i18rIr1C6dZ+1OtQLpTUj7W7+N2OXA5FL6AmJiE3kG2qqSGgBu11m8qpQYQadbzEnAFsFlrfYdS6lbgVuCWLI2x12GX3HbD+rdZdv6EdhPKzM8unV0eF39vaG3Z+9ctmEqBx2krRP0hI+naN214hzu+fRJup8PqFnXjmeP41uSR3POdk9Fa09AS4NqZY5kfU/kxIsS3c89FJ/PwFafidCg+q2vB0DpOQa144QPLQQ2RlYPX4ySYomnOiOJIhdZY01SqFUhXo4jsFPOKuRVW4pw4o4VcJ72FuDuI1nqv1vrN6OtG4H3gaOB84NHoYY8C38rG+HorbTl0E7clxtMbRqRh/bihRVbjeoiP9Q+GDVoCYdvYf6dSttc+qtjL8EH5bFg0jSevOY2vnjCU7z60jZkrtjD/16/RGjRo8odsPzt8UB75bifNgTDHlBbgcjq48pHXuXj1Nm7e8A6XTx8Tl1Pw0GWnMKQwL2U+xz/rWlgyaxxup8Pab4atxp5n1fzKLoX5GoZm3+FWW+W4aOZYSyGbJi9ByEWy3k9BKXUskUJ7fweGaa33QkRxKKWGpvjMAmABwOjRozMz0Bwl1lShVKRoXGKkTUsgXgGYM+raRr9l4mj0h+J6EZimmdjcA6dDMSDflZSL8OClkSY4dk5pQ4NDR3IsnE7FJQ/9PUlgPnzFqbafPdgcjHMAr5pXad3f9t0N3PXCDpadP4GxZYW4nA6cChp8AbTWPHb1FHYdaLF6QJv3c8/FJzO0KM/K99i+u4FHX/mUNVdNscJWh3TBxGOuEJpTKDizWqskuAm5TlaVglKqCNgI/JvW+rBSHfshaq1XA6shkqfQcyPMbexMFSvnRdpnm5E2pk8hNrpnzVVT2H/Yb2v/31Pvs5ytphN22MB8nrtuOr6AQX1zgKOKvay9+jRChsG+Q60UF7i5b/NHcTb/s8uHct0Z47ji4dfilEdZQn0hMwEuMcnsR+eemGRSWri2iie+dxrVexvZU++jtsnP8EH5hAzNJb961bYg34OXTqbJH7LMNwV5TlwuB+OHDWD9wml80RDp+3zj+rfZvrvB6jPQWWLNb21Va5UENyHXyVqPZqWUG9gEvKC1/p/oth3AzOgqYQSwRWs9vq3z9LfktVhSJZetXziNkKFxKuKij8yoIq21lUUc+7nYhjQAf7lpJl6PixKv24rWSVxJ1Db5uePbJ3HrM+/y5DWn8VFNs2Xzv/KR15OuEZuTEHvdlVs+ZsmscYwa7GX3QR/HDing9LuS+zCbNYvMSCqnA867P/IMVs2vjIskij3/sk3VSfb8dCagmc1P7HpEiE9ByFFyp0dztEHPr4H3TYUQ5XngcuCO6N/nsjC8rNPRrOVUPgStNaMHx3dhiu0BcPeFE9s0cYBpYnJRnO9if2OynfyWje+w7PwJDC704FCwYm4FdU0Bjh1SwK4DLRR47Md27JD4nIRYgTmkyINDqYj/IEUfZrfLEecg/7y+xTomVSTRicMjSWaJzzGdCWimH8M0a5mtPY8u9pLndnD/JZMk+kjoFWTLfDQDmA+8q5R6K7rt34kog/VKqauBfwIXZml8WaMjs1dTaYS15uErTuX37+5lVvkwir1uWgLhpAJuEB+ZlCrixvQ9mPb7QXlOPqhpwh+0F/DHlRVy+2+rqW0McPM547nuySPd0lbG2P9jr7H/sJ+nFkylviXI4Gif5J9/52QONvsBrMY5Z5cP5cF5lVwbU0F1xdwKmlpDDCnU1rOITbBLdV9ejytlklm6EtBi61KZWdkPXXYKwwbmR8Za2O1LCEJGyJr5KF30NfNRe/10E5WGabuPTQwzk8liZ6SxvV3tTBw/v/hknI5IW0unQ1HocREyNBetejUpwcsc09MLp+F0KoIhw9Yc9fj3TuPSX/09yeR090UT+aLBx5pXd/Gjc09kV10Lxw4pYP6v4xuSn10+lKWzv8wXDT6rV3Ntkz+u6U2qhj/dNQd1hVzpiicIHSR3zEe9mZ7+4bfXMyExF2FO5ShLIZjHXvPYG0mCE45kCMdG7hxTWoChId+tCIQ0d/z+fctJ/djVU9hTn7ra6G3Pv8cNZ42P66kQO2bAyito8AUtH8RndS2MLi1gTuUo6poD/P7dvVx7+tikc7xYXcOt557Ixau3xW2PbXqTaAJyOx2sXzCVkIZ8tyNlCY6eQMpeCH0BUQqdIBOVMdvrp5uoNFLZ0U3BaY7Z0Jr7L5lEfXOQAo+TlkCYEcV51Db6ufHpt+OEfW1jgO27G9h1oMXWTj58UD4/eGI723c3UL23kacXTrMd895DrQwpyrMqlsauFpbMGmdVUD33pBHWtRLPYSQsZO2a3rTVN3lIDrbHFIRcJivJa72VVO0w05mMZNqmE5OyzGSqxMSs2OSySaOKWTW/kg2LpuFUR1Y116x5g8O+IMGQwdLn3uPi1dtY+tx7BEKaP72/j6Wzy1m3YCpLZ5fz6CufsmjmWADu3byTlfMqLcWwbFM1eW4nz1btYfvuBiaNKuaOb59EWBs8cOnkuDEvn1PB8t9/wH1/+pCnFkzl2cXTWTq73OrF8Pt39zJ8UD4jS7wcO6SAezfvTEoiWzmvkrARjttm1/Smre+mwRegttHP5/Ut1Db6rVWTIAj2yEqhE2SiHWZ7ETGJjXY2Vu3mgUsnc/+fdsbVLDKVyeCCyEoibOiknsf3/2mnbYXQgfmR/xa1TX5GDMpj3YKphAyNy6F46R97ufuPO5P8EmeXD+Xx753G4dYQXzT4LOEPcP2ZxzMoGtn0k/PKGeh141AKrcETnfbXNvmt1YjpMB8xKI9DvpBlgkrV9CbVd1NWlMfB5oDVk6Elmhl9bGmh2PoFIQWiFDpBe6adrpDKR9FWtEycDd3lwB8Mc/M5J1qJYnBkprw+atpxOx1JQnNO5ShLIZifuWXjkQzjhy47hZKCvLiop9PGljGy5DMWzRwbVyPpxeoaqvc2suaqKVYjG4g4i0NhzWWPRaqlzpt2jOVQNstkv7mz1kp8M5v3rJhbwSFfiNElBQzId7frw7H7bv796ydS2+i3CviZ5y0ucDNYzEqCYItEH3WCdPsU0nU+w9DsaWjhK3dusbaZHcxOGD4Aj1PRGjKoawpQ1xywhPaGRdOYu/LVpPO9/MOv4nE6bBvomEqsJRDiqyu2JH12w6JpFBd4uOnpt6lt8vPYVVOszOSXbviKbULbw1ecyooXPuDmc06krslvRRqVDfDwk/MmoLVu16lv9yyfumYq33koOSpq3YKpHF1SYHseQehHSPRRd+loslNHI5Tsun7d89KOdrt+2Y3L63ZZM+VY045dmObyORU8+sqnDI4pf2EyssTLrgMtTDh6kO2YzVVMbaN94xzTv/Jf53854uxWRyKRnA77onlOh+LF6hqu/pfjrEijSaOKuXz6GC5a9WqHFKbdd+ML2tchCvfueZAg9CjiaO4kplA8uqSAsgHJ4Y7mjPWCB7YyY/nLXPDA1kj0j42D0zAMLp8+hmWbqrl49TaWbarm8uljbLt+GYZu02Ea66C++ZzxBEIGd184kTvnVvDw1k+TTEQ3fe0Eq9hdoiP3mNKCdquElhZ6WDW/Msm5vLFqN3XNAa59/E38obAVVQQQNrRt9VJzu0MpVs2vZN2Cqdw5N6K4OuPUT/xu8t321VLz3fLfXhBSISuFNJMqCiaxfwFAWBNnlzcFdmLXr46YmcyZ8nPXTWdvg5+bNsSHgZphpuZ1VLSi6LCB+Txz7XRaQ2EcSuF2KPI9DvYdbiUYNnA7HQwtysPlihekDofi6OJ81lw1hYPNEbPUo698yuXTx3DXCzvYU+9j+KB8bnr6HatZzkN/+SSpA9sDl05mwxv/5JErT6U1aLBsfXWb4+6MU39IYV6cU17CVAWhfUQppJnORChprVPWLoqlM4qmuTVsdSszjzUrnprF7kaWRIrOXX/m8YwuKcDhUJbSMZ3BsYJ75bxKThg2AJfLEWcaC2sdKWR35jgg4rg2o47MHgu1TX58gfCRHIeBeZGCfWEDp0OR53Lwva98CY3miodfaXfcsU799sx06axtJAj9BVEKaaYzEUodPbajiqauOUBNo9/2WNMclFjh9JnF0wmHNW6ng8eumkK+x8ltCe02F62tYt2CqUDEBPTT31bzYnUND19xKueeNIK9Da2EDCOuH8OKuRUcaAqwYm4FP/vd+9Zs33T0ej2uOAEdW9gu1bhj8zU66qSXLGNB6BxiXE0z7SWfdeXYVJ3E7JRHXXPA9tgRg/KtBDVzNr+n3keLP8yFq17lzP/5M/N/8xq7DjSz+PQvWZ3XICKc9x5qZcbyl7nkV3/n8uljmDSqmHs372R0aQE/+937eD1Olp0/gXULprLs/AmUDchj9GAvd/5hR1zLzD31kX4Nib6WVPc4fFA+W285nWcXz4gT+JlIJBSE/ogohTQTa7KwE2ZdObYzymNj1W7b9pJej5Mbn36bhY9Vxc3aPz3QHCdYb9rwDvXNQe6cW2EphtioItOss2jmWLbvbuCwL8CSWePQOnKc2fpzcKEHh8NBbZM/boxmw5lEIW53j6vmVzJ8QL6tUz8TiYSC0B8R81Gaac/Obbe/PfOGw6EYV1Zk2eJdUeevnfK44azxVpiree5BXhdFHneS03XlvEqW/u97cefYUx/J/j3kC/LDr42Pcx7HHlPsdTNpVDGhsE5KDjumtMDKOk68pmm6Ms9jCvHO2v97IpFQEARJXksr7dm529oPxCmLEm+kW5pZ+fNwa5ArHn693Zj9UMigtslPa8hg14Fmq0fxQ5edwriyIuucSik+r/dxw/q3kgTrsvMnEAhHfARPXjOVZZv+kdQXYenscjxOh6UQYvc9s3g6Qwrz4npHKwUf7G1M6v0wcdSgLmUXZ6I4YVeREtpCL8H2P6UohTTSXi+EVPufWTyduqYjNvKF/3os3zx5pBVFZM7AY+3zsec1MQXlvkOttsLa7D5mZiR/0eDD7XTwb+veiruO1+Pk/z5fzfbdDTy7eDqtwXBc8tuq+ZUMKfTgDxtxWdQmf//RGRxsCSYJ7KI8p2Wiit0eq6w6I0RzUfjmsrIShAQko7mnac/OnWp/a9CwhMikUcVcPOUYqwOZecxNG+LDM1NFH12z5g3u++4k61yLZo61sqUNw0gSWPdfMomfX3wyQ4ry0Gj2HWq1FMLIEi81jX5WbvnYij6KFb61jX7OLh/KnMpR1jU2Vu0mrLF1Aq9fOC2pKN89L+3g+jOPZ+FjVZ0WorkYWdSZ8GFByEXE0ZxG2osSSrXfGVMKYtHMsRxsDtgqj8QeynbRR9OPK2VIUR6/XfIv/OS8L8dlSx9oCmBozdLZ5UwaVcyeeh/XPbGduuYAP/tdNc3+yIrAVAjL51RY3c48LmeSw7fE62bJrOPjrrFk1vFx9xM7/mDYsC3KZyoE87jeHEUkDnChtyNKoYvYlZ1oL0oo1X6v54iyKPa6U4aVxvVQnl+ZFH3k9TiZN+0YvvvQNnyBcFIF1IVRc9SyTdX88GvjLcUwbmgRPzlvAicMG8Azi6fz55tm8siVUxiY72LJrHGsuWqKbUhtvS+YlCi3aG0V4WgkUuL4XQ6VtL200NOnhGhHw4cFIVcR81EXSGU3HldWxLCBkf4DYZt2kG1FEZlROqYJJrH95QOXTqbZH2LDommUDcjjqIH5SY7p1qBBIGRY5qBUq43YbOFlm6rZWdPEhKMG4nJFxrujKfne7Eg1K9Za20YdrXnlU6vkhbl96IC8PhVFlNjvoq08FUHIRcTR3AVSOYyf+N5pXBLTqD7RNt6WE9IwNDVNfkJhAw08sW0Xk48tpbTQQ2mhh8OtQVqDBg2+IJWji6lNcExfMvVYahv91DUH2Fi1m1vPPTHOL2GOMdYvsWHRNPwhg0df+dSqzNqes7wjz8F0aO873MoXDb64ct1nlw+NK4dd4nWzs7apTzlmc9EBLgg2iKM5XaSaIceWmLBzMKZyQj5/3Qz2H/bHCcZV8ysZXOjhg72N3P7b+DIR6xdOi3NMf73iaC6NUUbL51Twlx37efDSyrj+yA/Oq+SxV3ZZ5xnkdbPihQ+44azx1ky2MzbxtmbFDodCa53Ur+HF6hpu+6aO62fQ1+oT5aIDXBA6iiiFLpAqcSrROZooTFMJXF8gnKQsFj5WxTOLpzN8UL6VFWwK3dhCeotmjk3ZPW3FCx/E9Wq4b/OHzKkcxSuf1LFqXiUD813cfkFFnBDuTFJYewlnHT2XCFFByB1EKXQBuxnyqvmV/OKPH8YdlygAUwnJcIpqqcGQESd0vR4nIUPTGgzz8BWncu/mnZaPIPGzZuOa2KQzgP/4RjlPL5xKWVF+Ujls897WXDWFz+pa4voap7KJtyXQxb4uCL0PUQp03gZsN0Mu8bq54azxVO9tTCkAUwlJsxmM3YzaFLp2/ogVcysIho0Un3XYbv+ktpnhg/IZNjB1tziXQ/Hka5/xYnWNNcau2MmldLUg9D76vaM5nRmoHRGcdscA1hjKivJYMmscY4YUUpDntKKXTKduWVGelZDWEghz0tED2d/oj0v+inREy6f6i8akNpxmyexYX4fdMzCPNZ3DiQlmq+ZXMn7oANvVhiAIvQIpc2FHZ6JtehLD0DT4AuxtaLXyCWIV1N5DkUQzs/eyJZznVTLQ6+KjmuY4c8/okgL+Wd+CUlBz2E+DL2hFAAFsveV0y9mb6hmYkUqr5ldavRJi9z/xvdMYGW3SIwhCr0Oij+zIlQxUh0MRNrAUgjkOM4LJ43KyZNa4pPadC9dWsez8CVz5yOvWuUyldmxpIfsOt3Lj02+36exN9QzMDOpUCWY1jX68HlevchJLuKggtE2/X/unykB1uxxJGcs9TVsKqrTQw5ghhbb7zR4GiZ9xOBTDB+a324sh1TNo8AXjEswS95vCNZeJzTyvaWxlV10zFzywlRnLX05q9CMIgiiFlKUnmlpDGRcebZVIcDgUBXn2+83yF4mfgY418knV4ObkkYN4dvEMjhpk+imO7F8+p4KNVbtzOvPY9JWY3+O3H3iF/YdbKSuKrGx6e50lQegJ+r1PAZJNCk4HnHd/5v0MXe3HkOdyWNnLXXWUt2dWCYUMvjgUMRmZWdM3nDU+pzOP2/OVmMT6VwShH5E7PgWl1G+A2UCN1npCdNtgYB1wLLALuEhrXZ+J8STG2qdqIt/TppL2QjhT7Qe6HfbZXgKZy+VgZEkBXo+LEYPymTy6Iuft8e35SqB311kShJ4gW47mR4D7gTUx224FNmut71BK3Rp9f0sWxtblVo8ddWK2dVx7wjnV/kw4ezuSeZxLjtxU32NstVlJphOEeLJmPlJKHQtsilkp7ABmaq33KqVGAFu01uPbO09PFMTrTNtMU+h1NN+hL3fmyrV7SzWeYQPz8AWyr7QEIcvkVp6CjVJo0FoXx+yv11qXpPjsAmABwOjRoys/++yztI+vvSSzRKFX1xzoUL5DruRF9AS5eG+5tHIRhBzD9ofQK6OPtNartdanaK1PKSsr65FrmKaS2G5jqaqcmkKnI36IXMiLsGsQlA5y4d4SsfseBUFITS4lr+1XSo2IMR/VtPuJDNOW0Eu0X08aVcySWeMI64gANmeoXfVXQHp8Fj1p4unOvQmCkBvk0krheeDy6OvLgeeyOBZb2sojiI31nzSqmJvPGc/S597jK3duictzaK9lZyoSY+5T5U60d5y52ikrymPV/EruvnAi+w610GnSyAAACzFJREFU0uDrfqx+V+9NEITcISs+BaXUk8BMYAiwH7gN+F9gPTAa+Cdwodb6YHvn6oqjuat25vZi9WPPe/HqbSlt6125fkft9e0d93l9S8oaSieOGNjt1YLY8AWh15A7eQpa6++m2DWrp6/dVfOJYeiktpGr5lXypSGFSeGk7eU5dKWpjGm6mjSq2KqS2uALYhiG7XGprt1WDaV0OISlYY4g9G5yyXyUEdpyFnf2cwvXVrH3cGuSCactM1NX8bicnF0+lB9+bTzLNlVz8eptLNtUzYHmQNz127t2WzWUcr2OkSAIPU+/UwpdjZBpqy9zokLpCdt6aaGH//xGefIM/7GquOu3d+22aiiJQ1gQhFyKPsoIXY2Qaasv84hB+UnHDxuYx7oFUwlryHc7rGY5XcXhUDgdql2F1pFuZ0MK86RNpiAItvQ7pdDVvsGlhR5Wza+M6z62fE4Fj77yKZNHV1jHpfJZDCnsvp29owqtI6UypE2mIAh29MsqqT0VfQRdz+rtaCvPXCojIQhCryZ3oo96Kx2pFNoVn0VHhb3M8AVB6Gn6nVLo7my7PdNMZ30WhqHZd7jVNiLKbnUhIZ+CIPQk/S76qKshqR0lMbP54StOZe3Vp6HRKbOPv2jw9WiIaE/VOhIEoe/R71YKPV20zTTxPH/dDPY2tLJwbVXKFYmpoJbOLrddXbhd3dfZ4ocQBKEz9LuVQjoSy9qbeTscirCBpRDAfkViKqiVWz5m+ZyKuNyCFXMrcKVBaPf0ykgQhL5Fv1spdDUk1aSjM++OrEhMBbV9dwN3vbCDpbPLKS30MMjr5uYN73D/JZOgsHv3m4vlrAVByF363UohNoJn6y2n8+ziGZ0ypXR05t2RFUms/2H77gaWbarGHzK4ecM71Db505Jh3BMlNwRB6Lv0u5UCdC+Cp72Zd2y+wRPfO42f/raaF6trbFckpoJ6ZvF0WvxhPj3QzF0v7KC2yZ+2DOPurowEQehf9Eul0B3aCjm1My2tml/JsvMn4HA4bHMKHA7F0AH5GIWawjwX918yKa35B5LbIAhCZ+h35qPu0lbBOdtKqo9V4XA42m0F2ZNtI6UlpSAIHUVWCp2krZl3b3HqSiMcQRBSIUqhC6TySfSGHsWStyAIQluI+ShNGIZGo1l79Wk8fMWpTBpVnJNOXclbEAShLWSlkAZsHczzKhlenIdCsfeQL2fMNL3FxCUIQnaQlUIaSNWqs7k1zHn3b2XG8pe54IGt7NjfmPW6Q5K3IAhCW4hSSANtterMNTNNT7QKFQSh7yDmozTQVqvOWHLBTCN5C4IgtIWsFNKA3ex71fxKNlbtjjsuV8w0krcgCEIqZKWQBuxm3yVeNzecNZ7qvY1SXkIQhF5Dv+zRnCkkSUwQhBxGejRnmv7YOlMUoSD0bkQpCGlDsqUFofcjjmYhbUi2tCD0fmSl0A/IlElHsqUFofcjSqGPk0mTTm8oCCgIQtuI+aiPk0mTjmRLC0LvJ+dWCkqpc4BfAE7gV1rrO7I8pF5NJk06ki0tCL2fnFIKSikn8EvgLGAP8LpS6nmtdXV2R9Z7ybRJpz+G4QpCXyLXzEdTgI+01p9orQPAU8D5WR5Tr0ZMOoIgdIacWikARwOxBYP2AKdlaSx9AjHpCILQGXJNKdhJqqQ6HEqpBcACgNGjR/f0mHo9YtIRBKGj5Jr5aA8wKub9SOCLxIO01qu11qdorU8pKyvL2OAEQRD6OrmmFF4HximlxiilPMB3gOezPCZBEIR+Q06Zj7TWIaXUdcALREJSf6O1/keWhyUIgtBvyCmlAKC1/h3wu2yPQxAEoT+Sa+YjQRAEIYuIUhAEQRAsRCkIgiAIFqIUBEEQBAtRCoIgCIKFKAVBEATBQpSCIAiCYCFKQRAEQbAQpSAIgiBY5FxGcy5jGJq65oCUoBYEoc8iSqGDGIZmx/5Gq9+x2axm/LABohgEQegziPmog9Q1ByyFAJE+x9eseYO65kCWRyYIgpA+RCl0kEAoHNfnGCKKIRAKZ2lEgiAI6UeUQgfxuJxWn2OTkSVePC5nlkYkCIKQfkQpdJDSQg8PXXaKpRhMn0JpoSfLIxMEQUgf4mjuIA6HYvywATy7eIZEHwmC0GcRpdAJHA5F2YC8bA9DEAShxxDzkSAIgmAhSkEQBEGwEKUgCIIgWIhSEARBECxEKQiCIAgWSmud7TF0C6VULfBZtsfRTYYAB7I9iBxCnscR5FnEI8/jCN19Fge01uckbuz1SqEvoJR6Q2t9SrbHkSvI8ziCPIt45HkcoaeehZiPBEEQBAtRCoIgCIKFKIXcYHW2B5BjyPM4gjyLeOR5HKFHnoX4FARBEAQLWSkIgiAIFqIUBEEQBAtRChlGKfUbpVSNUuq9mG2DlVIvKaV2Rv+WZHOMmUIpNUop9bJS6n2l1D+UUtdHt/fX55GvlHpNKfV29Hn83+j2MUqpv0efxzqlVL9p4qGUciqltiulNkXf9+dnsUsp9a5S6i2l1BvRbWn/rYhSyDyPAIkJI7cCm7XW44DN0ff9gRBwo9b6RGAq8H2lVDn993n4gTO01hOBk4FzlFJTgeXAPdHnUQ9cncUxZprrgfdj3vfnZwFwutb65Jj8hLT/VkQpZBit9V+Agwmbzwcejb5+FPhWRgeVJbTWe7XWb0ZfNxL58R9N/30eWmvdFH3rjv7TwBnAhuj2fvM8lFIjgW8Av4q+V/TTZ9EGaf+tiFLIDYZprfdCRFACQ7M8noyjlDoWmAT8nX78PKLmkreAGuAl4GOgQWsdih6yh4ji7A/8HLgZMKLvS+m/zwIiE4QXlVJVSqkF0W1p/61I5zUh6yilioCNwL9prQ9HJoT9E611GDhZKVUMPAucaHdYZkeVeZRSs4EarXWVUmqmudnm0D7/LGKYobX+Qik1FHhJKfVBT1xEVgq5wX6l1AiA6N+aLI8nYyil3EQUwuNa62eim/vt8zDRWjcAW4j4WoqVUuYEbiTwRbbGlUFmAOcppXYBTxExG/2c/vksANBafxH9W0NkwjCFHvitiFLIDZ4HLo++vhx4LotjyRhRG/Gvgfe11v8Ts6u/Po+y6AoBpZQXOJOIn+VlYG70sH7xPLTWP9Jaj9RaHwt8B/iT1vpS+uGzAFBKFSqlBpivgbOB9+iB34pkNGcYpdSTwEwiZW/3A7cB/wusB0YD/wQu1FonOqP7HEqpfwH+CrzLEbvxvxPxK/TH51FBxFnoJDJhW6+1/i+l1HFEZsuDge3APK21P3sjzSxR89EPtdaz++uziN73s9G3LuAJrfXtSqlS0vxbEaUgCIIgWIj5SBAEQbAQpSAIgiBYiFIQBEEQLEQpCIIgCBaiFARBEAQLUQqCIAiChSgFQRAEwUKUgiCkiWg9+1dttp+klAoqpS7JxrgEoTOIUhCE9PE3YLJSKs/cEC3l8QDwitb6iayNTBA6iFRJFYT0sRXwECkBvi267TIiRe0mZ2tQgtAZZKUgCOljGxAmogSIFre7E7hfa/1uNgcmCB1FlIIgpIlo17S3iSoF4HYihf5uy9qgBKGTiFIQhPSyFZiqlJoMLAJu0lofzvKYBKHDSJVUQUgjSqmLgHXAP4CDWuuvZHlIgtApxNEsCOlla/TvCYhzWeiFiFIQhPTSBASAB7XW72R7MILQWcR8JAhpRCl1N3AJcILW+lC2xyMInUVWCoLQTZRSBcBE4F+B64m0RBSFIPRKRCkIQvc5k0jD9M+B67XWz7ZzvCDkLGI+EgRBECwkT0EQBEGwEKUgCIIgWIhSEARBECxEKQiCIAgWohQEQRAEC1EKgiAIgoUoBUEQBMHi/wNP/euIHWlaSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sns.scatterplot(model.y, model.y_hat)\n",
    "ax.set_xlabel(r'$y$', size = 16)\n",
    "ax.set_ylabel(r'$\\hat{y}$', rotation = 0, size = 16, labelpad = 15)\n",
    "ax.set_title(r'$y$ vs. $\\hat{y}$', size = 20, pad = 10)\n",
    "sns.despine()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
